{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TensorFlow version: 2.0.0\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "import datetime, os\n",
    "#hide tf logs\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'  # or any {'0', '1', '2'},\n",
    "#0 (default) shows all, 1 to filter out INFO logs, 2 to additionally filter out WARNING logs, and 3 to additionally filter out ERROR logs\n",
    "import scipy.optimize\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker\n",
    "import time\n",
    "from pyDOE import lhs         #Latin Hypercube Sampling\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import codecs, json\n",
    "\n",
    "# generates same random numbers each time\n",
    "# np.random.seed(1234)\n",
    "# tf.random.set_seed(1234)\n",
    "\n",
    "print(\"TensorFlow version: {}\".format(tf.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Data Prep*\n",
    "\n",
    "Training and Testing data is prepared from the solution file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1 = np.linspace(-1,1,256)  # 256 points between -1 and 1 [256x1]\n",
    "x_2 = np.linspace(1,-1,256)  # 256 points between 1 and -1 [256x1]\n",
    "\n",
    "X, Y = np.meshgrid(x_1,x_2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Data\n",
    "\n",
    "We prepare the test data to compare against the solution produced by the PINN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_u_test = np.hstack((X.flatten(order='F')[:,None], Y.flatten(order='F')[:,None]))\n",
    "\n",
    "# Domain bounds\n",
    "lb = np.array([-1, -1]) #lower bound\n",
    "ub = np.array([1, 1])  #upper bound\n",
    "\n",
    "a_1 = 1 \n",
    "a_2 = 4\n",
    "\n",
    "usol = np.sin(a_1 * np.pi * X) * np.sin(a_2 * np.pi * Y) #solution chosen for convinience  \n",
    "\n",
    "u = usol.flatten('F')[:,None] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trainingdata(N_u,N_f):\n",
    "    \n",
    "    leftedge_x = np.hstack((X[:,0][:,None], Y[:,0][:,None]))\n",
    "    leftedge_u = usol[:,0][:,None]\n",
    "    \n",
    "    rightedge_x = np.hstack((X[:,-1][:,None], Y[:,-1][:,None]))\n",
    "    rightedge_u = usol[:,-1][:,None]\n",
    "    \n",
    "    topedge_x = np.hstack((X[0,:][:,None], Y[0,:][:,None]))\n",
    "    topedge_u = usol[0,:][:,None]\n",
    "    \n",
    "    bottomedge_x = np.hstack((X[-1,:][:,None], Y[-1,:][:,None]))\n",
    "    bottomedge_u = usol[-1,:][:,None]\n",
    "    \n",
    "    all_X_u_train = np.vstack([leftedge_x, rightedge_x, bottomedge_x, topedge_x])\n",
    "    all_u_train = np.vstack([leftedge_u, rightedge_u, bottomedge_u, topedge_u])  \n",
    "     \n",
    "    #choose random N_u points for training\n",
    "    idx = np.random.choice(all_X_u_train.shape[0], N_u, replace=False) \n",
    "    \n",
    "    X_u_train = all_X_u_train[idx[0:N_u], :] #choose indices from  set 'idx' (x,t)\n",
    "    u_train = all_u_train[idx[0:N_u],:]      #choose corresponding u\n",
    "    \n",
    "    '''Collocation Points'''\n",
    "\n",
    "    # Latin Hypercube sampling for collocation points \n",
    "    # N_f sets of tuples(x,t)\n",
    "    X_f = lb + (ub-lb)*lhs(2,N_f) \n",
    "    X_f_train = np.vstack((X_f, X_u_train)) # append training points to collocation points \n",
    "    \n",
    "    return X_f_train, X_u_train, u_train \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PINN \n",
    "\n",
    "$W \\in \\mathcal{R}^{n_{l-1}\\times{n_l}}$ \n",
    "\n",
    "Creating sequential layers using the $\\textit{class}$ tf.Module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sequentialmodel(tf.Module): \n",
    "    def __init__(self, layers, name=None):\n",
    "\n",
    "        self.W = []  #Weights and biases\n",
    "        self.parameters = 0 #total number of parameters\n",
    "\n",
    "        for i in range(len(layers)-1):\n",
    "\n",
    "            input_dim = layers[i]\n",
    "            output_dim = layers[i+1]\n",
    "\n",
    "            #Xavier standard deviation \n",
    "            std_dv = np.sqrt((2.0/(input_dim + output_dim)))\n",
    "\n",
    "            #weights = normal distribution * Xavier standard deviation + 0\n",
    "            w = tf.random.normal([input_dim, output_dim], dtype = 'float64') * std_dv\n",
    "\n",
    "            w = tf.Variable(w, trainable=True, name = 'w' + str(i+1))\n",
    "\n",
    "            b = tf.Variable(tf.cast(tf.zeros([output_dim]), dtype = 'float64'), trainable = True, name = 'b' + str(i+1))\n",
    "\n",
    "            self.W.append(w)\n",
    "            self.W.append(b)\n",
    "\n",
    "            self.parameters +=  input_dim * output_dim + output_dim\n",
    "            \n",
    "        self.X = np.zeros(self.parameters) #store iterates\n",
    "        self.G = np.zeros(self.parameters) #store gradients\n",
    "        self.store = np.zeros((max_iter,2)) #store computed values for plotting\n",
    "        self.iter_counter = 0 # iteration counter for optimizer\n",
    "    \n",
    "    def evaluate(self,x):\n",
    "        \n",
    "        #preprocessing input \n",
    "        x = (x - lb)/(ub - lb) #feature scaling\n",
    "        \n",
    "        a = x\n",
    "        \n",
    "        for i in range(len(layers)-2):\n",
    "            \n",
    "            z = tf.add(tf.matmul(a, self.W[2*i]), self.W[2*i+1])\n",
    "            a = tf.nn.tanh(z)\n",
    "            \n",
    "        a = tf.add(tf.matmul(a, self.W[-2]), self.W[-1]) # For regression, no activation to last layer\n",
    "        \n",
    "        return a\n",
    "    \n",
    "    def get_weights(self):\n",
    "\n",
    "        parameters_1d = []  # [.... W_i,b_i.....  ] 1d array\n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "            \n",
    "            w_1d = tf.reshape(self.W[2*i],[-1])   #flatten weights \n",
    "            b_1d = tf.reshape(self.W[2*i+1],[-1]) #flatten biases\n",
    "            \n",
    "            parameters_1d = tf.concat([parameters_1d, w_1d], 0) #concat weights \n",
    "            parameters_1d = tf.concat([parameters_1d, b_1d], 0) #concat biases\n",
    "        \n",
    "        return parameters_1d\n",
    "        \n",
    "    def set_weights(self,parameters):\n",
    "                \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            shape_w = tf.shape(self.W[2*i]).numpy() # shape of the weight tensor\n",
    "            size_w = tf.size(self.W[2*i]).numpy() #size of the weight tensor \n",
    "            \n",
    "            shape_b = tf.shape(self.W[2*i+1]).numpy() # shape of the bias tensor\n",
    "            size_b = tf.size(self.W[2*i+1]).numpy() #size of the bias tensor \n",
    "                        \n",
    "            pick_w = parameters[0:size_w] #pick the weights \n",
    "            self.W[2*i].assign(tf.reshape(pick_w,shape_w)) # assign  \n",
    "            parameters = np.delete(parameters,np.arange(size_w),0) #delete \n",
    "            \n",
    "            pick_b = parameters[0:size_b] #pick the biases \n",
    "            self.W[2*i+1].assign(tf.reshape(pick_b,shape_b)) # assign \n",
    "            parameters = np.delete(parameters,np.arange(size_b),0) #delete \n",
    "\n",
    "            \n",
    "    def loss_BC(self,x,y):\n",
    "\n",
    "        loss_u = tf.reduce_mean(tf.square(y-self.evaluate(x)))\n",
    "        return loss_u\n",
    "\n",
    "    def loss_PDE(self, x_to_train_f):\n",
    "    \n",
    "        g = tf.Variable(x_to_train_f, dtype = 'float64', trainable = False)\n",
    "\n",
    "        k = 1    \n",
    "\n",
    "        x_1_f = g[:,0:1]\n",
    "        x_2_f = g[:,1:2]\n",
    "\n",
    "        with tf.GradientTape(persistent=True) as tape:\n",
    "\n",
    "            tape.watch(x_1_f)\n",
    "            tape.watch(x_2_f)\n",
    "\n",
    "            g = tf.stack([x_1_f[:,0], x_2_f[:,0]], axis=1)\n",
    "\n",
    "            u = self.evaluate(g)\n",
    "            u_x_1 = tape.gradient(u,x_1_f)\n",
    "            u_x_2 = tape.gradient(u,x_2_f)\n",
    "\n",
    "        u_xx_1 = tape.gradient(u_x_1,x_1_f)\n",
    "        u_xx_2 = tape.gradient(u_x_2,x_2_f)\n",
    "\n",
    "        del tape\n",
    "\n",
    "        q = -( (a_1*np.pi)**2 + (a_2*np.pi)**2 - k**2 ) * np.sin(a_1*np.pi*x_1_f) * np.sin(a_2*np.pi*x_2_f)\n",
    "\n",
    "        f = u_xx_1 + u_xx_2 + k**2 * u - q #residual\n",
    "\n",
    "        loss_f = tf.reduce_mean(tf.square(f))\n",
    "\n",
    "        return loss_f, f\n",
    "    \n",
    "    def loss(self,x,y,g):\n",
    "\n",
    "        loss_u = self.loss_BC(x,y)\n",
    "        loss_f, f = self.loss_PDE(g)\n",
    "\n",
    "        loss = loss_u + loss_f\n",
    "\n",
    "        return loss, loss_u, loss_f \n",
    "    \n",
    "    def optimizerfunc(self,parameters):\n",
    "        \n",
    "        self.set_weights(parameters)\n",
    "       \n",
    "        with tf.GradientTape() as tape:\n",
    "            tape.watch(self.trainable_variables)\n",
    "            \n",
    "            loss_val, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "            \n",
    "        grads = tape.gradient(loss_val,self.trainable_variables)\n",
    "                \n",
    "        del tape\n",
    "        \n",
    "        grads_1d = [ ] #store 1d grads \n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            grads_w_1d = tf.reshape(grads[2*i],[-1]) #flatten weights \n",
    "            grads_b_1d = tf.reshape(grads[2*i+1],[-1]) #flatten biases\n",
    "\n",
    "            grads_1d = tf.concat([grads_1d, grads_w_1d], 0) #concat grad_weights \n",
    "            grads_1d = tf.concat([grads_1d, grads_b_1d], 0) #concat grad_biases\n",
    "        \n",
    "        return loss_val.numpy(), grads_1d.numpy()\n",
    "    \n",
    "    def optimizer_callback(self,parameters):\n",
    "                \n",
    "        loss_value, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "        \n",
    "        u_pred = self.evaluate(X_u_test)\n",
    "        error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)\n",
    "        \n",
    "        tf.print(loss_value, loss_u, loss_f, error_vec)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(50):\n",
    "    \n",
    "    N_u = 400 #Total number of data points for 'u'\n",
    "    N_f = 10000 #Total number of collocation points \n",
    "\n",
    "    # Training data\n",
    "    X_f_train, X_u_train, u_train = trainingdata(N_u,N_f)\n",
    "\n",
    "    layers = np.array([2, 50, 50, 50, 1]) #3 hidden layers\n",
    "\n",
    "    maxcor = 200 \n",
    "    max_iter = 5000\n",
    "\n",
    "    PINN = Sequentialmodel(layers)\n",
    "\n",
    "    #load saved weights \n",
    "    new_checkpoint = tf.train.Checkpoint(model = PINN)   # get checkpoint of new model\n",
    "    chkp_path = \"Non_Stiff_problem/my_checkpoint\"                # define stored path\n",
    "    new_checkpoint.restore(chkp_path)                    # overwrite new model variables \n",
    "\n",
    "    init_params = PINN.get_weights().numpy()\n",
    "\n",
    "    start_time = time.time() \n",
    "\n",
    "    # train the model with Scipy L-BFGS optimizer\n",
    "    results = scipy.optimize.minimize(fun = PINN.optimizerfunc, \n",
    "                                      x0 = init_params, \n",
    "                                      args=(), \n",
    "                                      method='L-BFGS-B', \n",
    "                                      jac= True,        # If jac is True, fun is assumed to return the gradient along with the objective function\n",
    "                                      callback = PINN.optimizer_callback, \n",
    "                                      options = {'disp': None,\n",
    "                                                'maxcor': maxcor, \n",
    "                                                'ftol': 1 * np.finfo(float).eps,  #The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol\n",
    "                                                'gtol': 5e-10, \n",
    "                                                'maxfun':  50000, \n",
    "                                                'maxiter': max_iter,\n",
    "                                                'iprint': -1,   #print update every 50 iterations\n",
    "                                                'maxls': 50})\n",
    "\n",
    "    elapsed = time.time() - start_time                \n",
    "    print('Training time: %.2f' % (elapsed))\n",
    "    \n",
    "    print(PINN.loss(X_u_train, u_train, X_f_train))\n",
    "\n",
    "    PINN.set_weights(results.x)\n",
    "\n",
    "    ''' Model Accuracy ''' \n",
    "    u_pred = PINN.evaluate(X_u_test)\n",
    "\n",
    "    error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)        # Relative L2 Norm of the error (Vector)\n",
    "    print('Test Error: %.5f'  % (error_vec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05928775510204082\n",
      "0.06225121630695361\n",
      "0.298452\n",
      "1.3981489922579782\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAayklEQVR4nO3de3Bc5Znn8e/T3ZIl28I2kjAGm5FDwISLbRIFMsAQBkxhLuGSCQlkTMKsCTNby82bhMtOJYFUbYUEJgmpTQBDMkkK1gkwZJNhZggJAwskMWBjY2OsYJab5QsWwthgW7akfvaP0y1162LJ6vecVsu/T1WXuo9Odz/ukn969Z73Yu6OiIhUnlS5CxARkZFRgIuIVCgFuIhIhVKAi4hUKAW4iEiFyiT5Zg0NDd7U1JTkW4qIVLzly5e/4+6NfY8nGuBNTU0sW7YsybcUEal4ZvbmQMfVhSIiUqEU4CIiFUoBLiJSoRLtAxcR6ezspLW1lY6OjnKXMurU1NQwffp0qqqqhnW+AlxEEtXa2kpdXR1NTU2YWbnLGTXcnfb2dlpbW5k5c+awnqMuFBFJVEdHB/X19QrvPsyM+vr6ffrLRAEuIolTeA9sXz8XBbiISIVSgAd27g+e5rN3/6ncZYjIfkAXMQNbs3F7uUsQkf2EWuAisl969NFHmTt3LnPnzuXEE08km82Wu6R9NmQL3Mx+ApwHbHH3Y3PHDgR+CTQBbwCfdfet8ZUpImPRLf+6hpcD/9V69CEH8I1PHTPkeVdffTVPPfUU06ZN6zm2detWpkyZErSeOA2nBf5TYH6fYzcCj7v7EcDjucciIhXjnHPOYfbs2Vx33XU9xxYtWtTvvJHsG/ytb32Lr3zlK7S0tAx6Toj9iIdsgbv7U2bW1OfwBcBpufs/A54Ebii5GhHZrwynpRyHP/7xj7g7mzZtIpOJYvDRRx+lpaWF2267jcsuu4yLLrqICy+8kC984QvcfffdbN26lfr6ehYuXMiCBQs4//zzWbp0Kbfffjvf/OY3mTRpEvPnz6euro4lS5bwpS99iZ///Ofs3LmTbDbLD37wAzZv3lz0uoWt/5EYaR/4VHfflLu/GZg62IlmdqWZLTOzZW1tbSN8OxGRcB588EGOPPJIMpkM7s727dtpaGhgwYIFfPWrX2XlypVceuml3HDDDWSzWbq6upg8eTJ/+MMfePHFF/n85z/PokWLyGQytLS0UF1dzTXXXMO8efOYNWsWn/zkJ0mlUuzatYvJkyezbds2gKLXLTW8IcAoFHd3Mxv0bwF3XwwsBmhubi79bwYRkRJdeumlLFy4kMWLF1NbW8uPfvQjVq1axZw5c4AoaC+44AIAvva1r3HHHXfQ1tbG+vXrefHFF7nooouAaOLNmWeeyYwZM7jqqqu48847WbduHXPmzGHp0qX88Ic/ZNy4cT3vW/i6IYw0wN82s2nuvsnMpgFbglUkIhKzE044gdWrVxcd27BhA/feey8NDQ2sW7eOWbNmAXDMMcdw++23097ezvHHH8/KlSs58sgjeeeddzj44IO54YYb6O7u5rDDDuOggw7iwQcf5JRTTmHq1KlcfvnlzJgxg9NPP5358+cXvW4INpyO9Fwf+CMFo1BuA9rd/VYzuxE40N2vH+p1mpubfazvyNN0478B8Mat55a5EpHRae3atXzkIx8pdxmj1kCfj5ktd/fmvucO2QduZkuAPwGzzKzVzBYCtwJnmtk6YF7usYiIJGg4o1AuHeRbZwSupeKFGBYkIjJcmokZUHdWAS4iyVGAB9Rd0AJXa1xE4qYAD6gws9UaF5G4KcADyhYkeLda4CISMwV4QIWNbuW3iMRNAR5QYbeJulBEJG4K8IAKL1xm1QQXkZgpwAMqbHRX4NrwIvuV/WJDBxk+XcQU2Uf/cSNsXj30efvi4OPg7KEnhw+0oUOlUQs8oKy6UEQqxkAbOlxxxRXlK2gE1AIPyIu6UBTgIkMaRks5DgNt6LBr1y7Wrl3LzTffTEtLC0uWLOGWW27p2cjh61//OuvXr+/ZvOGss87iiSeeKNqwIWlqgQdUOPJE+S0yeg20ocOKFSu4+OKLufnmm5k0aRIbN24s2sgBKNq84ZVXXum3YUPS1AIPSF0oIpVhoA0dnnvuOWbPng3Azp07+23kABRt3lBfX89dd91VtGFD0hTgARV1oSjARUatgTZ0uOuuu9i8eTMPPPAACxcuZPny5UUbOQBFmzecddZZ/TZsSJoCPKBs0WJWZSxERPbZPffcU/T4tNNO63fOt7/97aLHn/rUp+IsaUjqAw8oqxa4iCRIAR5QcR94GQsRkf2CAjwgTaUXkSQpwAMqXo1QAS4yGP3/GNi+fi4K8ICKN3QoXx0io1lNTQ3t7e0K8T7cnfb2dmpqaob9HI1CCUjjwEWGNn36dFpbW2lrayt3KaNOTU0N06dPH/b5CvCAFOAiQ6uqqmLmzJnlLmNMUBdKQO4D3xcRiYMCPCDNxBSRJCnAA9I4cBFJkgI8oKINHZTgIhIzBXhAxZmtABeReCnAg1IXiogkRwEeUPGmxkpwEYmXAjygrHbkEZEEKcADKsxsTRMWkbiVFOBmtsjM1pjZS2a2xMyGP4l/DCra0KGMdYjI/mHEAW5mhwLXAM3ufiyQBi4JVVgl0kQeEUlSqV0oGaDWzDLAeGBj6SVVLk3kEZEkjTjA3X0DcDvwFrAJ2Obuj/U9z8yuNLNlZrZsrK8+pha4iCSplC6UKcAFwEzgEGCCmS3oe567L3b3ZndvbmxsHHmlFaB4U2MFuIjEq5QulHnA6+7e5u6dwMPASWHKqkyFkZ3Vhg4iErNSAvwt4BNmNt7MDDgDWBumrMqkPTFFJEml9IE/CzwEvACszr3W4kB1VaTCVrcuYopI3ErakcfdvwF8I1AtFc/38khEJDTNxAxIwwhFJEkK8IA0jFBEkqQAD8i1oYOIJEgBHpAiW0SSpAAPKKthhCKSIAV4QEV94JrIIyIxU4AHpOVkRSRJCvCYqAtFROKmAA9Ii1mJSJIU4AEVjwMvXx0isn9QgAeU1UQeEUmQAjwgL+pCKWMhIrJfUIAHVBja6gMXkbgpwANytJiViCRHAR6Q+sBFJEkK8ICKu1DKV4eI7B8U4AFpLRQRSZICPKDCyFZ+i0jcFOABaVNjEUmSAjwgzcQUkSQpwANSH7iIJEkBHpAyW0SSpAAPqDC/s+pDEZGYKcADKr6IWcZCRGS/oAAPyDUTU0QSpAAPSBs6iEiSFOAB+SD3RUTioAAPSMMIRSRJCvCANJFHRJKkAI9BytQCF5H4KcADyo/9TqdMk3pEJHYlBbiZTTazh8ysxczWmtlfhiqsEuUzO2WmUSgiErtMic+/A3jU3T9jZtXA+AA1Vax8t0k6ZeoDF5HYjTjAzWwScCpwOYC77wH2hCmrMuUb3Wkz9YGLSOxK6UKZCbQB/2xmK8zsXjObEKiuipTvNkmpD1xEElBKgGeAjwJ3uvvxwA7gxr4nmdmVZrbMzJa1tbWV8HajnwNmGoUiIskoJcBbgVZ3fzb3+CGiQC/i7ovdvdndmxsbG0t4u9Ev607KDFMXiogkYMQB7u6bgfVmNit36Azg5SBVVSh3MKIWuPJbROJW6iiUq4H7cyNQXgP+rvSSKpeDWuAikpiSAtzdVwLNYUqpfNlcE1wtcBFJgmZihuRReKfUAheRBCjAA8q6Y1guwMtdjYiMdQrwgDzXAjcNIxSRBCjAA8o6mBmmPnARSYACPCDHcxN5tJiViMRPAR5Q7zhw9YGLSPwU4AG5e08XivrARSRuCvCAook8+S6UclcjImOdAjygbK4FrsWsRCQJCvCAXBN5RCRBCvCAoguX+bVQyl2NiIx1CvCgPJrIAxpGKCKxU4AHlM3mNnRIaSKPiMRPAR6QE23ooD5wEUmCAjygbG4ij5nRrfwWkZgpwAPy3Foo0XrgSnARiZcCPKBoJqYm8ohIMhTgARXuSt+tcYQiEjMFeEDuuogpIslRgAeULViNUPktInFTgAeU35U+ldJaKCISPwV4QL270qsLRUTipwAPyaPw1looIpIEBXhA0a70aBy4iCRCAR6Q51rgKTO6FeAiEjMFeEDZgok82Wy5qxGRsU4BHlC2YCq9LmKKSNwU4EG59sQUkcQowAOKWuBoHLiIJEIBHlB+Kr3pIqaIJEABHpCm0otIkhTgAUWrERppXcQUkQSUHOBmljazFWb2SIiCKlnheuBaTlZE4haiBX4tsDbA61S8nok8KXWhiEj8SgpwM5sOnAvcG6acylY4lV4tcBGJW6kt8O8D1wODzjs0syvNbJmZLWtrayvx7Ua3fAs8ndJqhCISvxEHuJmdB2xx9+V7O8/dF7t7s7s3NzY2jvTtKkJ+OVmtRigiSSilBX4ycL6ZvQH8AjjdzO4LUlWFijZ0QFPpRSQRIw5wd7/J3ae7exNwCfCf7r4gWGUVyN0xjLQ2dBCRBGgceECem0pvGkYoIgnIhHgRd38SeDLEa1WybG4qfVrDCEUkAWqBBxTNxNQwQhFJhgI8oGzBRB4tZiUicVOABxStRphfzEoBLiLxUoAH1NMHrouYIpIABXhA2Wx+Q4doIo9a4SISJwV4QL3LyRqAZmOKSKwU4AHl+8DTuU9Vk3lEJE4K8ICyBVuqgYYSiki8FOAB5Tc1TqfyXSgKcBGJjwI8oGhHnt4+cLXARSROCvCACnfkgWhUiohIXBTgAWXzFzGj/NZsTBGJlQI8oGy0n0NPH7i6UEQkTgrwgBwv7kJRC1xEYqQADyiaiamLmCKSDAV4QL0TeRTgIhI/BXhAfceBK8BFJE4K8IAKd+QBjUIRkXgpwAPKerQSoVrgIpIEBXhAPX3guogpIglQgAfUrwtFAS4iMVKAB5TfE1MBLiJJUIAHlHUvGoXSpQAXkRgpwAPKL2aVSUUfq2ZiikicFOABZftM5OnqVoCLSHwU4AHlL2Jm0uoDF5H4KcADimZi9l7E7NSC4CISIwV4QPlx4Jn8KBR1oYhIjBTgAXVnvegipkahiEicFOAB5afS5/vAu9SFIiIxUoAH4rkhg1pOVkSSMuIAN7MZZvaEmb1sZmvM7NqQhVWafFanzKjKdaF0qg9cRGKUKeG5XcCX3f0FM6sDlpvZ79z95UC1VZT8pB2D3i6UbnWhiEh8RtwCd/dN7v5C7v77wFrg0FCFVZp8d0kqZVSl8y1wBbiIxCdIH7iZNQHHA88O8L0rzWyZmS1ra2sL8XajUn7WfDplVOVa4OpCEZE4lRzgZjYR+BfgOnff3vf77r7Y3ZvdvbmxsbHUtxu1ugsuYmbS+WGEaoGLSHxKCnAzqyIK7/vd/eEwJVWmni4UUwtcRJJRyigUA34MrHX374YrqTLlhxGmU0Z1rgW+u0stcBGJTykt8JOBy4DTzWxl7nZOoLoqTmEL3CwKcV3EFJE4jXgYobs/QzRqTigYB56bxFOVNvaoBS4iMdJMzECyBRcxAaozKQW4iMRKAR5IvgslvyP9uEya3V3d5SxJRMY4BXggPS3wXBO8pipFR6da4CISHwV4IPkh3ym1wEUkIQrwQLp7hhFGj2uqUuxSC1xEYqQAD6RwGCFAbXWajj1qgYtIfBTggWQLJvIAjK/OsLOzq5wlicgYpwAPpCs3bT7TE+Bpdu5WC1xE4qMAD6R3HHgU4HU1Gd7frRa4iMRHAR5Ivg88v5lDXU0V73d0lrMkERnjFOCBdPW5iDmptoqOziwdnepGEZF4KMAD6XsRc/L4KgDe26lWuIjEQwEeSM9U+lyA10+oBqB9x+6y1SQiY5sCPJC+a6HUTxwHQPsHe8pWk4iMbQrwQLr6tMAbcwG+5X21wEUkHgrwQLpzi6Hk98NsqIsCvE0BLiIxUYAH0nciz4TqNLVVad75QAEuIvFQgAfS1WccuJnRUFetABeR2CjAA+kJ8FTvR9owcZwCXERiowAPpCu3gXG+CwWgfsI4jUIRkdgowAPp24UC0DCxmnd3KMBFJB4K8EB6L2L2fqQHTogC3HOzNEVEQsqUu4CKt/NdeO0Jjl37DHdUtVD/y+9B106onsDnd6Y5PJ2h4w9vUjvzBJh6HGSqy12xiIwRCvCR+GALrLwfWv4NNiwHz3KsZZhkB2JVs6DuIOjcwfgd73Jq6k1qf/909LzqiXDUeXDcxfCh0yCtj19ERk4JMlzu8MYz8Py90PIIZLvgkI/CqdfDh+dx17o6vvPY/+PPl80nk0kDsPqVNr74k2f59YK/YI69Cq8+Dmt/A6t+ARMOgo8vhOaFMLGxzP84EalECvChZLOw7rfw9D9B6/NQOwVO/Af42OXQcETPabtbXgGgOt3bBx4taGVsopE5xxwHx1wE5/4TvPp7WP4zePJb8PR3YfZn4ZRFUH94wv84EalkCvDBZLPQ8q/wf78Db78Ekw+Dc78Lc/8Wqmr6nb67K0t1OoVZ4SiUaDp90UiUzDg46tzo1vYKPHsnrPzfUZfMMZ+Gv/oyTD069n+eiFQ+BXhf+eB+8tuwZQ3UfxguuhuO/cxe+6w7OrupqSoe1HNgbknZQddDaTwSzvsenHYT/Ol/wfM/hpceivrJ/+q/w6EfC/bPEpGxRwGe190Jqx+CZ74H7/w5Cu5P3wPH/g2k0kM+PQrw4vOqMymmjK9iy/sdPcd+vXIDHZ3dfO7jh/WeOPEgOPObcPJ18Ozd8OxdUT/7h/4aTroaDj8dClr2IiKgAI+GAa64D56/B957Cw46Bv7mx1F/9TCCu+dl9nRTW93//KkH1PD29qgF3tHZzbW/WAnA2cdN44CaquKTxx8If30TnHRV1Bpf+iO479PQeBSc+PfR6JVxdSP+p4rI2LJ/Bng2C+uXwotLYNWD0LULDjsJzv4OHDl/RK3dD3Z3UVfT/+M8eFINm7fvAuDlTdt7jq9u3cbJH24Y+MXG1cEp18En/ius+VUU5I8sgt/+I3zkfJhzCcw8dZ9+wYjI2LP/BHjXHmh9Dl75Lbz0MGxvhUwtzL4YTrgSDj6upJd/v6OTieP6f5zTp9Sy4q33cHfWbOwN8DUb9xLgeZlxUVjP/lw0Ambl/fDSr6JhiLUHwhFnwpFnQdOpGooosh8qKcDNbD5wB5AG7nX3W4NUFcIHbbBpJWxcGYXfG89A5w5IZeDwM2DezTDrbBg3McjbtX+wh6MPOaDf8Q81TGTbrk7aPtjNS63bmDy+itqqNKs3bB/gVQZhBjNOiG7zb41+Cb3yKKx7DFb9Mjpnyszo+9PmQMMsaJwFBxwKKa2WIDJWjTjAzSwN/BA4E2gFnjez37j7y6GK67FnJ3TuhK4O6OyIgnj3+9CxHTrei2ZG7miD7Rth6+vw7mvQsa33+fVHwNxLo4uBTadAzaSg5XV2Z2l9bxdnHj213/eOmx6913Ovv8sfX3uHjx02hdrqNEtfa6erO9uzg8+wVdXCMRdGt2w3bHgB3voTrH8WXnuyN9AB0tVQNw0mTYeJU2F8fXSrnRzNCh03MfqaqYleNzMuek66GtJVkKqKfuGl0tHNCr5aKvrFoourMhD33C3b59Ydfc129x7rud9dcN97z+255V6PwdYWyv88Wu7nM3dLpQd4nO7/vZ77hedb79dRqJQW+AnAq+7+GoCZ/QK4AAgf4A9cFk1+2Zuq8VB3MExpiob81R8etUYPng01/VvGIb3+zg4MmDtjcr/vHT9jModOruXqJStwh+vPOooJ49I8smoTVy9ZwZ0LShgqmErDjI9Ht7wd7dEomrYW2PombN8A2zbA5lWwsx12bR35+w3GUhT/58l9heL7PY97HuzlNUfnf5gxadDF1nyAc3yA+wN8HZMG+Bnv+Tntc7/nKQX3//6posl/QSoa6Up5ZvYZYL67X5F7fBlwortf1ee8K4Ercw9nAX/ex7dqAN4ZUZHloXrjVWn1QuXVrHrjNZJ6/8Ld+13oiv0iprsvBhaP9PlmtszdmwOWFCvVG69Kqxcqr2bVG6+Q9ZZyhWsDMKPg8fTcMRERSUApAf48cISZzTSzauAS4DdhyhIRkaGMuAvF3bvM7Crgt0TDCH/i7muCVdZrxN0vZaJ641Vp9ULl1ax64xWs3hFfxBQRkfLSLA8RkQqlABcRqVAVEeBmdpuZtZjZKjP7lZlNLndNAzGz+Wb2ZzN71cxuLHc9e2NmM8zsCTN72czWmNm15a5pOMwsbWYrzOyRctcyFDObbGYP5X5215rZX5a7pr0xs0W5n4WXzGyJmfXfuaTMzOwnZrbFzF4qOHagmf3OzNblvk4pZ42FBqk3WJ5VRIADvwOOdffZwCvATWWup5+CpQXOBo4GLjWz0by1ThfwZXc/GvgE8N9Geb151wJry13EMN0BPOruRwFzGMV1m9mhwDVAs7sfSzQw4ZLyVjWgnwLz+xy7EXjc3Y8AHs89Hi1+Sv96g+VZRQS4uz/m7l25h0uJxpyPNj1LC7j7HiC/tMCo5O6b3P2F3P33icLl0PJWtXdmNh04F7i33LUMxcwmAacCPwZw9z3u/l5ZixpaBqg1swwwHthY5nr6cfengHf7HL4A+Fnu/s+AC5OsaW8GqjdknlVEgPfxX4D/KHcRAzgUWF/wuJVRHoh5ZtYEHA88W+ZShvJ94HogW+Y6hmMm0Ab8c67L514zm1Duogbj7huA24G3gE3ANnd/rLxVDdtUd9+Uu78Z6L+q3OhVUp6NmgA3s9/n+t763i4oOOcfif70v798lY4tZjYR+BfgOnffhzVuk2Vm5wFb3H15uWsZpgzwUeBOdz8e2MHo+tO+SK7f+AKiXzyHABPMbEF5q9p37pWzmlaIPBs1Gzq4+7y9fd/MLgfOA87w0Tl4veKWFjCzKqLwvt/dHy53PUM4GTjfzM4BaoADzOw+dx+tIdMKtLp7/q+ahxjFAQ7MA1539zYAM3sYOAm4r6xVDc/bZjbN3TeZ2TRgS7kLGkqoPBs1LfC9yW0ccT1wvrvvLHc9g6iopQXMzIj6Z9e6+3fLXc9Q3P0md5/u7k1En+1/juLwxt03A+vNbFbu0BnEsdRyOG8BnzCz8bmfjTMYxRdd+/gN8MXc/S8Cvy5jLUMKmWcVMRPTzF4FxgHtuUNL3f0fyljSgHKtw+/Tu7TA/yxvRYMzs1OAp4HV9PYp/w93//fyVTU8ZnYa8BV3P6/MpeyVmc0luuBaDbwG/J27x7AgexhmdgvwOaI/61cAV7j77vJWVczMlgCnES3J+jbwDeD/AA8AhwFvAp91974XOstikHpvIlCeVUSAi4hIfxXRhSIiIv0pwEVEKpQCXESkQinARUQqlAJcRKRCKcBFRCqUAlxEpEL9f1CW8F5osUajAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s = np.array(np.loadtxt('50_runs_stiff_transfer.txt'))\n",
    "t = np.array(np.loadtxt('50_runs_stiff.txt'))\n",
    "\n",
    "print(np.mean(s))\n",
    "print(np.std(s))\n",
    "\n",
    "sns.distplot(s, hist=False,\n",
    "                 kde_kws={\"shade\": False},\n",
    "                 norm_hist=False, label=r'$\\varepsilon_{transfer}$')\n",
    "\n",
    "print(np.mean(t))\n",
    "print(np.std(t))\n",
    "\n",
    "sns.distplot(t, hist=False,\n",
    "                 kde_kws={\"shade\": False},\n",
    "                 norm_hist=False, label=r'$\\varepsilon_{base}$')\n",
    "\n",
    "# plt.savefig('Transfer_histogram.png', dpi = 500, bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
